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Abstract 

In this paper we present a Maple library (MOPs) for computing Jack, Hermite, 
Laguerre, and Jacobi multivariate polynomials, as well as eigenvalue statistics for 
the Hermite, Laguerre, and Jacobi ensembles of Random Matrix theory. We also 
compute multivariate hypergeometric functions, and offer both symbolic and nu- 
merical evaluations for all these quantities. 

We prove that all algorithms are well-defined, analyze their complexity, and 
illustrate their performance in practice. Finally, we also present a few of the possible 
applications of this library. 
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1 Introduction 

1.1 Motivation 

There is no need for us to review the impact that classical orthogonal polynomial 
and special functions theory has had for applications in mathematics, science, en- 
gineering and computations. By the middle of the last century, handbooks had 
been compiled that could be found on nearly everyone's bookshelf. In our time, 
handbooks join forces with mathematical software and new applications making 
the subject as relevant today as it was over a century ago. 

We believe that the modern day extension of these scalar functions are the mul- 
tivariate orthogonal polynomials or MOPs along with their special function coun- 
terparts. 

The multivariate cases are far richer, yet at this time they are understudied, 
underapplied, and important applications may be being missed. Algorithms for 
their computation have not been studied systematically, and software suitable for 
scientific computing hardly exists. At this time there are no handbooks, no collection 
of software, and no collection of applications, though since April 2004 entries are 
being introduced into Eric Weinstein's Mathworld website 1 . 

Development of such software may thus be seen as a whole area of research 
ripe for study. This paper might be thought of as a first step in this direction; 
undoubtedly better software will emerge in time. 

We recall that scalar orthogonal polynomials are defined by a positive weight 
function w(x) defined on an interval I cK. We define the inner product 

<f,9>w = J f{x)g(x)w(x)dx 

and the sequence of polynomials Pq{x) , p\{x) , P2(x) , . . ., such that pk(x) has degree 
k, and such that < Pi,Pj > w = if i ^ j. This sequence is the sequence of orthogonal 
polynomials with respect to the weight function w(x). 

There is a (scalar) complex version of this inner product (/ C C) where we use 
g instead of <?; this induces a different set of orthogonal polynomials. 

We now define the multivariate version of the inner product, and the correspond- 
ing orthogonal polynomials. We take any weight function w(x) defined on a segment 
/, and create an n-dimensional weight function which is symmetric in each of its n 
coordinates, and incorporates a repulsion factor which depends on a "Boltzmann" 
constant f3 (or a temperature factor a = 2/f3) which is not seen in the univariate 
case: 

n 

W( Xl ,...,x n )= 11 Ix.-xj] 2 ^ Y[w( Xl ) . (1) 

l<i<j<n i—1 
Mathworld, URL http://mathworld.wolfram.com/ 
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We define multivariate orthogonal polynomials p™{x\, . . . ,x n ) with respect to 
the weight W(x\, . . . ,x n ). The polynomials are symmetric: they take the same 
value for any permutation of the n coordinates a^, and they satisfy 

/n 
p"{xi, . . .,x n )p™(x!, . . . ,x n )Y]_\xi - Xjf Y[ w(xi)dxi . . . dx n = 5 Kfl , 

i<3 3 = 1 

where k represents the "multivariate degrees" of p™ (the exponent of the leading 
term). 

We begin with our fourth example: symmetric multivariate Hermite polynomi- 
als. We take w(x) = e~ x / 2 , so that the integral is over all of K n . The polyno- 
mials are denoted H"(x). Our second and third examples are w(x) — x a e~ x and 
w(x) = x ai (l — x) a2 . These are the Laguerre L"' a and Jacobi J^ a ^- a ^ polyno- 
mials. Special cases of the Jacobi polynomials are the Chebyshev and Legendre 
polynomials. 

Our first example, the Jack polynomials, generalizes the monomial scalar func- 
tions, x k . These polynomials are orthogonal on the unit circle: w = 1 and I = the 
unit circle in the complex plane. Therefore I n may be thought of as an n dimen- 
sional torus. The orthogonality of the Jack polynomials may be found in formula 
(10.35) in Macdonald's book page 383]. 

Tables ^ |5J |3 and 0] give the coefficients of the Jack, Hermite, Laguerre, and 
Jacobi in terms of the monomial symmetric functions (for the first) and the Jack 
polynomials (for the last three). We take all degrees up to total degree 4 for the 
Jack polynomials, up to total degree 3 for the Hermite polynomials, and up to 
degree 2 for Laguerre and Jacobi; the coefficients can be seen by a simple call to 
the procedures, e.g. 2 , 

> jack(a, [2],' J'); 

> hermite(a, [1, 1, 1], n,' C'); 

> laguerre (a, [1, l],g,n, 'C'); 

> jacobi (a, [l],c/i,g 2 ,n, 'C'); 



1.2 History and connection to Random Matrix Theory 

The Jack polynomials have a very rich history. They represent a family of 
orthogonal polynomials dependent on a positive parameter a, and some of 
them are more famous than others. There are three values of a which have 
been studied independently, namely, a = 2,1,1/2. The Jack polynomials 
corresponding to a = 1 are better known as the Schur functions; the a = 2 
Jack polynomials are better known as the zonal polynomials, and the Jack 
polynomials corresponding to a = 1/2 are known as the quaternion zonal 
polynomials. 

In an attempt to evaluate the integral (J2J) in connection with the non- 
central Wishart distribution, James JHj discovered the zonal polynomials in 

2 Note the use of a for a and g for 7 in the calls. 
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Table 1: Coefficients of the Jack "J" polynomial expressed in monomial basis 
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Table 2: Coefficients of the Hermite polynomial expressed in Jack "C" polynomial basis 
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Table 3: Coefficients of the Laguerre polynomial expressed in Jack "C" polynomial 
basis 
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Table 4: Coefficients of the Jacobi polynomial expressed in Jack "C" polynomial basis 
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1960. 

(tr(AHBH T )) k (H T dH) = Y,c K Z K (A)Z K (B) . (2) 

Inspired by the work of James JHj and Hua |12j . in his own attempt to 
evaluate ©, Jack was lead to define the polynomials eventually associated 
with his name ^21- More explicitly, Jack orthogonalized the forgotten sym- 
metric functions 26, page 22], using the inner product given in Definition 
12. 9P . He studied the resulting one-parameter (a) class of polynomials (which 
now bear his name), and for a = 1 he proved were the Schur functions, while 
for a = 2 he conjectured to be the zonal polynomials (and proved his claim 
in a very special case). 

He consequently generalized the a parameter to any real non-zero number, 
and noted that for a = — 1 he obtained yet another special class of functions, 
which he called the "augmented" monomial symmetric functions. Later it 
was noted that the orthogonalizing inner product was positive definite only 
if a > 0. 

During the next decade, the study of Jack polynomials intensified; Mac- 
donald (2H1 page 387] points out that in 1974, H.O.Foulkes [HI raised the 
question of finding combinatorial interpretations for the Jack polynomials. 
This question was satisfactorily answered in 1997 by Knop and Sahi |21| . 

In the late '80s, the Jack polynomials were the subject of investigation 
in Macdonald's book j2^ and Stanley's paper [HJ; these two authors gen- 
eralized many of the known properties of the Schur functions and zonal 
polynomials to Jack polynomials. 

As mentioned, an important application of the Jack polynomials came in 
conjunction with random matrix theory and statistics of the 2/a-ensembles. 
Below we mention a few of the researchers who have made significant con- 
tributions in this area. 

James |16| was one of the first to make the connection between the zonal 
polynomials (a = 2 Jack polynomials) and the 1-ensembles, when he calcu- 

3 The authors would like to thank Plamen Koev for an in-detail explanation of this fact. 
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lated statistical averages of zonal polynomials over the 1-Laguerre ensemble 
(Wishart central and non-central distributions). 

At about the same time, Constantine and Muirhead provided a general- 
ization of the hyper geometric series, using the zonal polynomials, and studied 
the multivariate Laguerre polynomials for a = 1 (for a reference, see |28|). 

In a survey paper, James defined and described multivariate Laguerre, 
Hermite and Jacobi polynomials for a = 1 JH]- Chikuse [3| studied more 
extensively the multivariate Hermite polynomials for a = 1. 

In the early '90s, Kaneko [20] studied the general a binomial coefficients, 
and used them in connection with the study of hyper geometric series and 
multivariate Jacobi polynomials. He also studied Selberg-type integrals and 
established the connection with generalized Jacobi polynomials. A few years 
later, Okounkov and Olshanski [22] considered shifted Jack polynomials for 
all a, and proved that they were the same as the generalized binomial coef- 
ficients. 

Kadell ^H] was perhaps the first to consider averages of many valued Jack 
polynomials, with his study of the average of the Jack polynomial of parame- 
ter 1/k (with k an integer) over the corresponding 2k- Jacobi ensemble. Later 
it was noticed that constraining k to be an integer was unnecessary. 

Lasalle (221 ESI > considered all three types of general a multivariate 
polynomials, and among many other things computed generating functions 
for them. 

The last results that we mention here are those of Forrester and Baker 
[2], who studied in detail the multivariate, general a Hermite and Laguerre 
polynomials, in connection with the 2/a-Hermite and Laguerre ensembles 
(some of their work built on Lasalle |23l I25j ) . For a good reference on mul- 
tivariate generalizations of many of the univariate properties of the Hermite 
and Laguerre ensembles, see (§]. 

2 Multivariate Functions: Definitions, Properties, 
and Algorithms 

2.1 Partitions and Symmetric Polynomials 

Definition 2.1. A partition X is a finite, ordered, non-increasing sequence 
of positive integers Ai > A2 > A3 > . . . > A/. 

Throughout this paper, we will refer to / = /(A) as the length of A, and 
to k = |A| = Y?i=i as the sum of A. 

Remark 2.2. Naturally, one can remove the constraint "finite" from the 
definition of the partition, and replace it with "of finite sum", since one can 
always "pad" a partition with Os at the end; in this context I becomes the 
index of the smallest non-zero component of the partition A. 

We will work with two orderings of the partitions. The first one is the 
lexicographic one, denoted by <. 
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Definition 2.3. We say that A < k in lexicographical ordering if for the 
largest integer m such that A« = k, l for all i < m, we have X m < K m . If 
\ m < K m , we say that A < k. 

Remark 2.4. This is a total ordering of the partitions. 

The second ordering is the dominance ordering, sometimes also called the 
natural ordering. 

Definition 2.5. We say that A ■< k (or, equivalently, that k "dominates" 
X) if, given m = m&x{length(K) , length(X)} , 



3 3 





< E Ki ' V j < m , 


and 


i=l 


i=l 




m 


m 






= J2 K i ■ 




i=i 


i=l 





If one of the inequalities above is strict, we say that A -< K. 

Remark 2.6. Note that we compare two partitions only if they sum to the 
same integer. Also note that even with this constraint, ^ is only a partial 
ordering of the set of partitions of a given number: for example, [4, 1, 1] and 
[3, 3] are incomparable. 

The above summarizes what the user should know about partitions in 
order to use our library. 

Definition 2.7. A symmetric polynomial of m variables, 
polynomial which is invariant under every permutation of x±, . . . , x m . 

Remark 2.8. The symmetric polynomials form a vector space over R. 

Over the course of time, combinatorialists have defined a variety of homo- 
geneous bases for this vector space; each such basis is indexed by partitions 
(which correspond to the terms of highest order in lexicographical ordering of 
the polynomial). By homogeneity we mean that all terms of a polynomial in 
the basis have the same total degree (but this degree varies from polynomial 
to polynomial). 

Some of these homogeneous bases are displayed in the table below: 



Name 


Definition for I = 1 


Definition for I > 1 


power-sum functions 


P\i — Z^i 7=1 x j 


PX = UUlPXi 


elementary functions 


e Ai = 11 j : . j,. .... A x ji ■ ■ ■ x j Xl 


ex = n'=i exi 


complete homogeneous functions 


h\i = J2j 1 <j 2 <...<j M x h ■ ■ ■ x 3\ 1 


^ = nu hx z 



Another important basis is given by the monomial functions m, 

mx = E X l(l) X a(2) • • • X a(m) 5 
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here S\ is the set of permutations giving distinct terms in the sum; A is 
considered as infinite. 

The last basis we mentioned distinguishes itself from the other ones in 
two ways; the advantage is that it is very easy to visualize, and proving 
that it is indeed a basis is immediate. The disadvantage is that it is not 
multiplicative 4 . 

Monomials seem to be the basis of choice for most people working in 
statistics or engineering. Combinatorialists often prefer to express series in 
the power-sum basis, because of connections with character theory. 



2.2 Multivariate Gamma Function 

Before we proceed, we will need to define the multivariate Gamma function 
for arbitrary a; the real and complex versions are familiar from the literature, 
and the arbitrary a > case represents an immediate extension: 

r«(o) = ^(m-l)^) g r Lill] . (3) 

i=l ^ a ' 

Just as the univariate Gamma function generalizes to the multivariate 
one, the shifted factorial (Pochhammer symbol, rising factorial) 

v !k r(o) 

becomes the generalized shifted factorial. We call 

length(K) . . lengthen) „ / i-1 , \ 

ck= n - n w 

the generalized shifted factorial, or generalized Pochhammer symbol. 



2.3 Jack Polynomials (the Multivariate Monomials) 

Jack polynomials allow for several equivalent definitions (up to certain nor- 
malization constraints). In addition to the definition presented in the intro- 
duction (at the end of Section we present here two more (Definitions 
12.91 and !2.1()|) . Definition 12 . 91 arose in combinatorics, whereas Definition l2.1()l 
arose in statistics. We will mainly work with Definition 12. 1U1 

Definition 2.9. (following Macdonald }26$ ) The Jack polynomials are 
orthogonal with respect to the inner product defined below on power- sum func- 
tions 

(P\,Pn)» = a liX) z x 5x^, 

4 For x £ {p, e,h}, x\x^ = x p , where p can be obtained in an algorithmic fashion from A and fi 
(sometimes by mere concatenation and reordering). In general, m\m^ is not a monomial. 
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where z\ = Yl a i-i ai > a % being the number of occurrences of i in A. In 
i=i 

addition, 

p \ = m A + Yl u %v m v ■ 

There are two main normalizations of the Jack polynomials used in com- 
binatorics, the "J" normalization (which makes the coefficient of the lowest- 
order monomial, [1™], be exactly n!) and the "P" normalization (which is 
monic, and is given in Definition 12. 9|) , To convert between these normaliza- 
tions, see Tables El and IH1 In Table 03 I m = (1, 1, 1, . . . , 1), where the number 
of variables is m. 

We use the notation k h k for k a partition of k, and for Y17=i — 
1 -§(*-!))• 

Definition 2.10. (following Muirhead, [28]) The Jack polynomial C% is 
the only homogeneous polynomial eigenfunction of the following Laplace- 
Beltrami-type operator 

" d 2 . 2 ^ x? d 



d* = x 2 — |- _ 



dx- dxj ' 

with eigenvalue + k{m — 1), having highest-order term corresponding to 
K. In addition, 

C°(xi,x 2 , • • -,x m ) = (xi + x 2 + . ..x m ) k . 

k h fc, !(ft)<m 

Remark 2.11. T/ie "C" normalization for the Jack polynomial allows for 
defining scalar hyper geometric functions of multivariate (or matrix) argu- 
ment. These are useful for computing Selberg-type integrals and other quan- 
tities which appear in various fields, from the theory of random walks to 
multivariate statistics and quantum many-body problems. 

Remark 2.12. David M. Jackson [1^1 pointed out that the D* operator also 
appears in algebraic geometry, for example in the context of ramified covers. 

Definition 2.13. Given the diagram of a partition k (see Figure 1), define 
a K (s) (the "arm-length" ) as the number of squares to the right of s; l K (s) (the 
"leg-length") as the number of squares below s; /i* (s) = l K (s) + a(l + a K (s)) 
(the "upper hook length") and /i£(s) = l K (s) + 1 + aa K (s) (the "lower hook 
length"). 

Finally, a further definition is needed in order to present the conversion 
table. 

Definition 2.14. Let 

c(K,a) = YlK(s) , 



c'(a,n) = Yl h *(s) 



j K = c{a,n) c'(a,n) 



9 



Figure 1: The Arm-length and the Leg- length. 

K= (4,3,1,1) 

a K (s)= 2 

'k(s)= 1 



where h* K and h% have been defined above. 

To explain the conversions between "J" , "P" , and "C" , we recall the def- 
inition of the generalized Gamma function and generalized shifted factorial 
from Section POl 

We can now present Tables El and EI the entries have been filled out using 
James ^7], Forrester and Baker and Stanley pi] . 



Table 5: Values of the different normalizations of Jack polynomials of partition k and 
parameter a at I m . 



Normalization 


Basic Property 


Value at I m 


= (1,1,1 


...,1) 


C 


sums to {x\ + X2 + . . . + x n ) k 


CZ(Im) 


_ a 2k k\ f 


m \ 
a I 


J 


has trailing coefficient n\ 


■Wm) 






P 


is monic 




c(a,n) 


V / K, 



Table 6: Conversions between the three normalizations for the Jack polynomials; the 
a(V, W) entry above is defined as V£[x\, . . . , x m ) =a(V, W)W"(xx, . . . , x m ). 
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J 


P 


c 




a k fc! 


a k fc! 


3k 




J 


3k 
a k fc! 




c(k, a) 


p 


c'(k,o) 


1 




a k fc! 


c(k,«) 
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2.4 Algorithm used to compute the Jack Polynomials 

From the Laplace-Beltrami equation, one can find an expansion for the Jack 
polynomials of the type 

C"(xi,X2, ■ ■ .,x m ) = ^2c^ fl m x (x 1 ,x 2 , ■ ■ .,%) , 

A<K 

where A and k are both partitions of the same integer \k\, and the order 
imposed on partitions is the lexicographic one. The coefficients c° x depend 
on all three parameters; m\(xi,X2, ■ ■ ■ ,x m ) is the monomial function corre- 
sponding to A. 

Note that as a consequence of the above, if 1(k) > m, C"(xi, x 2 , ■ ■ ■ , x m ) = 
("there is no highest-order term"). 
Using the eigenfunction equation 

D*C: = (p a K + k(m-l))CZ , (5) 

where 

m „ 
p« = J>(**-l--(i-l)) 

i=i 

one can obtain a recurrence for c" > from which the Jack polynomials can be 
explicitly calculated. This recurrence is 

2 

a 



^ Pk ~ Px 



« E ((h + t)-{l 3 -t))c^, (6) 



where A = (h, . . . , k, . . . , lj, . . . ,l m ), n = (h, . . . ,k+t, ...,lj-t,.. .,l m ), and 
fj, has the property that, when properly reordered, it is between A (strictly) 
and k in lexicographic order. 

In fact we can do better, using two propositions found in Macdonald's 
book PHI (10.13), (10.15)]. Roughly the content of the two propositions is 
that the Jack polynomials, in "P" normalization, can be written as 

K = m K + E <a™a , 

A-<K 

with u" x > whenever k y A (the ordering imposed on partitions here is 
the dominance ordering). 

Thus it follows that the recurrence can be improved to 

2 

a 



£ Uli + Q-ilj-t))^, (7) 



where A = (h, . . . ,k, . . . . . .,l m ), P = (h, ■ ■ -,k+t, ...,lj-t,.. .,l m ), and 

fx has the property that, when properly reordered, it is between A (strictly) 
and k in domination order. 
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This recurrence, at first glance, seems to be enough to compute all coeffi- 
cients c° A , once c" K is found. However, one has to account for the possibility 
that p° = for some A different from k; what can one do in that case? 

Fortunately, this never happens. We first need the following well known 
Proposition. 

Proposition 2.15. The dominance ordering is a lattice on the set of parti- 
tions of a given number. In particular, between any partitions K and A such 
that k y X, there exists a "path" on this lattice, a = k >~ a 1 >~ . . . >~ a 1 = A, 
such that a t+1 differs from a 1 in the following way: there exists i\ < i% such 
that a l+1 and a 1 agree in all places but i\ and i2, (cr i+1 )i 1 = (cr 1 )^ — 1, and 
(o-^) i2 = (a%+l. 

Now we can prove that we never divide by in computing Recurrence |7| 

Lemma 2.16. If A -< k, then ^ p", for all a > 0. 

Proof. Let A -< k be two partitions, let m = m&x{length(K),length(\)} , 
and assume that there is some a > such that 

Since the two partitions sum to the same number, the above is equivalent to 

m „ m 

»=i i=i 

The right-hand side is non-negative (as an immediate consequence of the 
strict ordering). 

We show that the left-hand side is positive by induction. For that we will 
use Proposition 12.15] which shows that it is enough to prove that 

m 

$>?-A?>0 

i=l 

in the case when k and A differ only in two places, i% < i^. Note that if 
= Ajj + 1 and n i2 = \ i2 — 1, this implies that > n i2 + 2. Hence 

m 

kf - Aj 2 = k 2 n - Xl + k 2 %2 - A? = 2k h - 1 - 2k l2 - 1 > 2 > , 

i=l 

and we are done. □ 

Proposition 12.151 ensures thus that once c" K is determined, every other 
non-zero coefficient is uniquely determined. 

Finally, for c" K we use the following formula (deduced on the basis of 
Table Eland the fact that P" has highest-order coefficient 1): 

_ a k k\ 
c'( K ,a) ' 
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Remark 2.17. It is worth mentioning that, from the recurrence^ by letting 
a — ► oo, the coefficient c" A goes to faster than c" K , for any A / k. Thus, 
at a = oo, the Jack "P" polynomial (which is monic) is the symmetric 
monomial. This could also be seen from the weight functions, as at a = oo, 
the "interdependence" term Yl 

I ' a (see for example^) disappears 

l<i<j<n 

and the variables separate. 

2.5 Multivariate binomial coefficients 

Many algebraic quantities (and the identities they satisfy) can be extended 
from the univariate case to the multivariate case through Jack polynomials. 
One such example is the multivariate, or generalized, binomial coefficient. 

Definition 2.18. We define the multivariate (or generalized) binomial co- 
efficients (Z) as 

C%(xi + l,x 2 + l,...,x m + l) _ A /k\ C%(xi,x 2 , ■ ■ .,x m ) 

where a C k means that ai < k% for all i. 

The generalized binomial coefficients depend on a, but are independent 
of both the number of variables m and the normalization of the Jack poly- 
nomials (the latter independence is easily seen from the definition). 

The multivariate binomial coefficients generalize the univariate ones; some 
simple properties of the former are straightforward generalizations of prop- 
erties of the latter. For example, 

((0)) =1, 

((i)) = \ R \ ' 

□ =0 if a 2^, 

□ =8 K \i\K\ = \a\, 

(*) ^ if \K\ = H + 1, iff<7 = K (i) , 

where Kr^ = (k\, . . . , ki — 1, . . . , k m ). The above are true for all k and a, and 
a subject to the constraints. 

2.6 Algorithm used to compute the multivariate binomial co- 
efficients 

One can prove, using the eigenfunction equation (J5J) and the definition of the 
generalized binomial coefficients, that 
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where \a\ = s, \k\ = k, erW = (a± . . . , o~i + 1, ...,a m ). All generalized 
binomial coefficients can be found recursively, once one has a way to compute 
the so-called "contiguous" coefficients ( CT ) . 

To compute the contiguous coefficients, we use Proposition 2 from [20] . 
applied to k = crW, and simplified slightly: 

(») 



a 



where g° \ is 



Here 



\sea- / \s&a / 



h°(s), if s is not in the ith column of a , 
/i*(s), otherwise. 

(s), if s is not in the zth column of a , 



CT< ' |^ h^ } (s), otherwise. 

Knowing the contiguous coefficients allows for computing all the generalized 
binomial coefficients. 

Remark 2.19. The generalized binomial coefficients are independent of the 
number of variables. They are rational functions of a. 

2.7 Multivariate Orthogonal Polynomials 
2.7.1 Jacobi Polynomials 

These polynomials represent the Gram-Schmidt orthogonalization of the 
Jack polynomials C" with respect to the Jacobi weight function 

J{ ' r^ + ^ + i) f}^ J 

x Y[ \ x i ~ Xj\ 2 / a dxi . . . dx m , (11) 

i<j 

on the hypercube [0, l] m . For the purpose of well-definitedness we assume 

51,92 >-l. (12) 



Define 



_d_ 

dxi ' 



E 
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then the Jacobi polynomials are eigenfunctions of the following Laplace- 
Beltrami operator: 

D* + (gi +g 2 + 2)E -5*- ( 9l + l)e , (13) 

with eigenvalue p° + \ n\(gi + 92 + ~ 1) + 2). 

2.7.2 Algorithm used to compute the Jacobi Polynomials 

Using the fact that the Jacobi polynomials are eigenfunctions of the operator 
()13|) . one obtains that these polynomials can be written in Jack polynomial 
basis as 

Ta,gi,ga( T . T \ - (n, , m ~ 1 ^ rj a (T \ C -1 )*^ C£(xi, . . . ,x m ) 

J K ( Xl ,...,x m )-( gi + a +^CAI m )^ {gi + ^ + l)a c „ {Im) 

where the coefficients c" CT satisfy the recurrence 

2 L(o)( (T ) C " ctW ' 



(92 + 91 + i(m - 1) + 2)(k - s) + p% - p% \ i allowable 



with the previous notation for p° and . The question is again whether 
the denominator is always nonzero. 

Proposition 2.20. Under these assumptions, (g2 + Q\ + n( m — -0 + 2)(fc — 
s ) + Pk ~ Pa i s never 0. 

Proof. The proof is very similar with the corresponding proof of Section \'2A\ 
the two crucial facts here are that one needs one show it for the case k = crW , 
and that gi and g2 are both larger than —1 (due to (|12|l). □ 

Letting c" K = 1 for all k and a allows all the coefficients to be uniquely 
determined. 

2.7.3 Laguerre Polynomials 

The multivariate Laguerre polynomials are orthogonal with respect to the 
Laguerre weight function 

dpt(xi,...,x m ) = — — t e-^' YlxJYl\x i -x j \ 2/a dxi...dx m ,(15) 

r^(7 + — + i) Y 

on the interval [0, oo) m . Note that for the purpose of well-definitedness, we 
must have 7 > — 1. 

This weight function can be obtained from the Jacobi weight function Q1U|) 
of the previous subsection by substituting (gi+g2+^(fn— l)+2)~ 1 (xi, . . . , x m ) 
for (a?!, . . . ,x m ) and then taking the limit as 52 — * 00. The same limiting 
process applied to the Jacobi polynomials yields the Laguerre polynomials. 
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Under the transformation mentioned above, the Jacobi differential oper- 
ator becomes 

8*-E+(-y + l)e , (16) 

and the Laguerre polynomials are eigenfunctions of this operator with eigen- 
value \k\. 

2.7.4 Algorithm used to compute the Laguerre Polynomials 

The Laguerre polynomials have an explicit expansion in Jack polynomial 
basis, which depends on the generalized binomial coefficients: 

T a ,l(^ t \ — U a. m-i , -,\ r a (T \ ("^Q C"(xi, . . . ,x m ) 

k jxi, . . . , x m j - ( 7 + — + i) K u K [i m ) — | m srjY) C a (I m ) ' 

Note that the coefficient of C"(xi, . . . ,x m ) in L"' 7 (xi, . . . ,x m ) is (— l) k . 

2.7.5 Hermite Polynomials 

The multivariate Hermite polynomials are orthogonal with respect to the 
Hermite weight function 

(r(i _i_ l)) m 

dn%( Xl ,...,X m ) = 2- m ' 2 n rn(m-l)/a-m/2 Li_Lg^ x (1 7 ) 

x e -^ x ? /2 Y[\xi - Xj \ 2/a d Xl . . . dxll§) 

on M. m . 

This weight function can be obtained by taking (y+y/jxi, 7+^/7x2, • • • , 7+ 
^/ 7 yx m ) in (|15[) . and then letting 7 go to infinity; note that the only remaining 
parameter is a. 

Under this limiting process, the differential operator becomes 

5** — E , (19) 

where 



^ — ' rtrf rv — ' T; — 



— ' dxj dxj 

i 1 ijtj J 

The Hermite polynomials are eigenfunctions of this operator with eigen- 
value \k\. 

Remark 2.21. Similarly, 

lim 7 ~ fc/2 ^' 7 (7 + y/yx!, 7 + Jyx 2 , . . . ,7 + Jyx m ) = (-l) k H%( Xl , ...,x m ) 

7-^00 
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2.7.6 Algorithm used to compute the Hermite Polynomials 

Using the corresponding Hermite differential operator (|19|) . we obtain the 
following recurrence for the coefficients of the polynomial. Let 




In the above, i < j take on all admissible values, and we choose c" K = 

Alternatively, we can obtain the coefficients directly through the limiting 
process described in Remark 12.211 



k+s 
2 



C-{X)]H-{X) = 

^cr V-'rr 



M-3 



where 



F(r, a, m, k, a) 



k + s 
J- 2 



—J 



F(r, a, m, k, a) ,(22) 



(r + ±(m + a-l)) K 
( r + I( m + a _l)) CT 



We use the above formula to calculate a single coefficient, n, for reasons 
of smaller computational complexity (in computing integrals with respect to 
the Hermite weight function; see Section EH)- 

Note that if a 2 K or k ^ s (mod 2), then the above is 0. 

2.8 Hypergeometric functions 

The hypergeometric functions are perhaps the easiest to generalize from 
univariate to multivariate. For the multivariate versions, a good reference is 
Forrester's unpublished book 9 . 

Definition 2.22. We define the hypergeometric function p F" of parameters 
ai, . . . , a p , respectively bi, . . . , b q and of variables (x\, . . . , x m ) by 



pFg (<2>l, . . . , Op! &1, . . . , bq] X\, . . . , X m ) — ^ ^ 



k=0 Khk 



(ai) K . . . (a p ) K 
k\ {b 1 ) K ...{b q ) K 



C K [x\ , . . . , x m ) 



Note that this is a formal definition; p < q is needed in order for the 
hypergeometric series to converge everywhere, and when p = q + 1, there 
is a nontrivial convergence radius. When p > q + 2, the series converges 
everywhere except at 0, with one notable exception, made by the polynomial 
hypergeometrics, i.e. those for which some ctj is a negative integer, which 
forces the series to terminate after a finite number of terms. 
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This definition of a hyper geometric function assumes an argument (xi, . . . , x r 
M. m ; similarly one can extend the definition to hypergeometric functions of 
arguments in (xi, . . . , x m ; y\, . . . , y m ; . . .) S W m x W" 1 x ... by inserting an 
additional 

C%(yi, . . . y m )/C"(l, . . . , 1) for each extra vector in R m . 

Hypergeometric functions provide answers to many statistics and statistics- 
related questions; below are two examples. 



1. Krishnaiah and Chang |22j have proved in 1971 that the density of the 
smallest root of a real (a = 2) Wishart matrix with m variables and n 
degrees of freedom such that p = is an integer is proportional 
to 

p(x) = x? m e~ xm ' 2 2 F (-p, ^±1 -2/ m _!/x) . 

Note that the joint eigenvalue density of the matrix described above is 
given by d/U^ with a = 2 and 7 = p. 

In [7j we extend this to any a and any positive integer 7 = p. We obtain 
that for this case the density of the smallest eigenvalue is proportional 
to 

p(x) = x? m e~ xm ' 2 2 Fg(-p, - + 1; -21^/x) . (23) 

a 

2. The largest eigenvalue (Ji) distribution for a Wishart real matrix with 
m variables and n degrees of freedom (a = 2, 7 = n ~™ -1 ) can be 
expressed as 

T m \Um+l)\ (X\mn/2 11 1 

P ^ h <X ^ = r ri fZ ZvH (?) iF 1 (-n,-(n+m+l);--xI m ) . 
r m |_^(n + m + 1)J V2/ 2 2 2 

The above is a corollary of a stronger theorem proved by Constantine 
H], and it can also be found in Muirhead [23 page 421]. 
This result generalizes to any a and 7 (as noted in [7j) to 

r m [i(m-l) + l] /xW7+(m-l)/a+lj 



r m [7 + ^(^-i) + 2] V2 

. (m - 1) 2 „ x 1 \ 

x 1^17+^ - + a + -(m-l) + 2; - -xl m ) 

a a 2 ' 

non-central real (a = 1) Wishart matrix A = Z'Z, with Z a matrix 
of independent Gaussians with mean M and variance J n xE, and with 
matrix of noncentrality parameters Q = M T M , the moments of the 
determinant 



3 Software 

3.1 The model for MOPS 

We have initially chosen as a model for MOPS the Symmetric Functions 
(SF) package by John Stembridge 5 . Our library is compatible with SF, and 



5 The SF package can be found at the URL http://www.math.lsa.umich.edU/~jrs/maple.html#SF 
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our procedures m2p, p2m, and m2m are designed to complement the SF 
procedures torn and top (for a comparison, see Section 0}. Though in time 
our vision of MOPS has changed, we are grateful to John Stembridge for 
his ground-breaking work. 

3.2 System requirements and installation guide 

Our library was initially developed for Maple 7, and later for Maple 8. Ex- 
perience has shown that it is also compatible with Maple 9. 

We have developed MOPS on various Unix machines and one Windows 
XP machine; the same version of the library is compatible with both operat- 
ing systems. Below we provide an installation guide, which can also be found 
on the MOPS webpage, located at http://www.math.berkeley.edu/~dumitriu/mopspage.html 

For Unix users: download the file mops. zip into your home directory. 
Unzip the file using the command 

unzip MOPS. zip 

This should create a new MOPS directory in your home directory; the 
new MOPS directory contains a subdirectory named Help Files and 4 files, 
named maple.hdb, MOPS.ind, MOPS.lib, and MOPS.mws. The last file 
contains an archive of all the procedures we wrote for MOPS. 

Return to your home directory (e.g. '/home/usr/vis/joesmith/'), and 
create a .mapleinit file; write in 

newJibname := '/home/usr/vis/joesmith/MOPS' ; 
libname := libname, newJibname ; 

and then save the .mapleinit file. 

All that is left is to call the library (in Maple) using the standard Maple 
command for libraries, i.e. 

> with(MOPS) ; 
each time you need to use the library. 

For Windows users: users: place the downloaded file in your C:\ directory 
(or in a more appropriate place of your choosing) . 

Unzip the file using Winzip; this should create a new C:\MOPS direc- 
tory in your home directory; the MOPS directory contains a subdirectory 
entitled Help Files and 4 files, named maple.hdb, MOPS.ind, MOPS.lib, and 
MOPS.mws . The last file contains an archive of all the procedures we wrote 
for MOPs. 
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1. In the \Maple\bin folder, create a .mapleinit file. In it you should 
write 

newJibname := 'C:\\MOPS' ; 

libname := libname, newJibname ; 

and then save the file. 

You will need to call (in Maple) the library each time you need to use 
it, using the standard command 

> with(MOPS) ; 

2. For some unknown reasons, the instructions in variant 1 do not always 
work on a Windows XP machine. In that case, you will have to type 
in the pathway for the library, each time you will need to use it. Upon 
opening a Maple window, you will have to type in 

> newJibname := 'C:\\MOPS' ; 

> libname := libname, newJibname ; 

> with(MOPS); 

each time you need the library. 



For Mac users: the instructions are similar to the instructions for Win- 
dows. 



Regardless of the Operating System, we suggest you perform a check 
before you start using the library. For example, if you type in 

> jack(a, [3], 2, 'P') ; 

the answer should be 

roi , 3m[2,l] 

m[3] + TT^T 



3.3 Routines: name, syntax, and description 

We give here a list of the 30 routines, and a brief description of what they do; 
for more extensive mathematical explanation we refer to Section[21 Note that 
most of them are set up to do calculations both symbolically and numerically. 
We use the following notations: 

1. k, a for partitions, 

2. k, s for the corresponding partition sums, 

3. n, m for the number of variables, 
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4. a. for the Jack parameter, 

5. i,j for the location of a square in a partition, 

6. A p , B q for two lists of real parameters, one of length p and the other 
of length q, 

7. I for a computational limit, 

8. 7,51,92 for additional Laguerre or Jacobi parameters, 

9. [x] for either a number or a list of variables, 

10. r for a real parameter, 

11. N for the normalization, 

12. exp for an expression. 
Some parameters are optional. 

Three of the routines, namely, m2p, p2m, and m2m, are alternatives to 
the routines torn and top from the SF package to the conversion between the 
monomial and the power sum basis. For a comparison between our routines 
and their SF counterparts, see Section 0J 

Some of the routines are used as building blocks for others; for a tree 
description of the dependencies, see Figure El 
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Procedure 


Syntax 


Description 


arm 


armf/t i i) 


the arm-length of a partition at a square 


conjugate 


conjugate(K) 


the conjugate partition 


expH 


expH(a, exp, n) 


the expected value of an expression in terms of 
monomials with respect to the 2/ct-Hermite 
distribution with n variables 


expHjacks 


expHjacks(a, exp, n) 


the expected value of an expression in terms of Jack 
nolvnomials with resnect to the 2 /ry- Hermite 
distribution with n variables 


expJ 


expJ(a,exp,g 1 ,g 2 ,n) 


the expected value of an expression 

in tprms of monomials with rpsnprt to thp 

2/ a, gi, (72-Jacobi distribution with n variables 


expJjacks 


expjjacks(a, exp, g\ ,g 2 ,n) 


the expected value of an expression 

in terms of Ja.ck nolvnomials with resnect to the 

111 U V_/l .1. .1. ll > V-/1 r Clival V l_/ \J 1. y J. 1 v_/ 1111 l_l 1 ij VVlLfll 1 ljU k-' U U v-' Ullv 

2/o, gi, g 2 - Jacobi distribution with n variables 


expL 


expL(a, exp, 7, n) 


the expected value of an expression 

in terms of monomials with respect to the 

2/q, 7-Laguerre distribution with n variables 


expLjacks 


expLjacks(a, exp, 7, n) 


the expected value of an expression 

in terms of Ja.ck nolvnomials with resnect to the 

ill vj v_' 1 lllkj v_/l * r LliOll i_/ \_j± y 1 1 v_/ 1 1 1 1 ell v_> yviLfll iv_>lj i_/ v_> v> u U v-' U llv 

2/a, 7-Laguerre distribution on n variables 


gbinomial 


^binomial ( cv k rr\ 


the generalized binomial coefficient 


ghypergeom 


ghypergeom(a, A p , B q , [x],'N', I) 


the generalized hypergeometric function 
rorrpsnonrlinp' frj naramptpr lists J^,. pvalnatpH 

\^ Wl 1 V_ O J-/ VJllLllllg Li W ^JCAiL CllllV LjV^I llu UO J. A£j . J_J Q . V V Cll LLC*j UVjU. 

at the point x £ W 1 (or written 
symbolically for n variables), or (optional) 
as a series "cut" after terms that sum up to I 


gsfact 


gsfact(a, r, k) 


the generalized shifted factorial (or generalized 
Pochhammer symbol) 


hermite 


hermite(a, k, [x],'N') 


the multivariate Hermite polynomial 
written in Jack polynomial basis 


hermite2 


hermite2fa k \x\ 'N') 


alternative wav of pomnutiner the multivariate 

C-ii yj v^i 1 1 i_i u i y \-y y y t_i y v^i Vjwiii i_/ i_i uiii t-. uiiv iimi ui v hi l < i Uv 

Hermite polynomial written in Jack polynomial basis 


issubpar 


issubpar (a, k) 


checks if a is a subpartition of k 


iack 


jack(a, k, [x],' N') 


the Jack polynomial as a linear combination 
of monomials 


jack 2 jack 


jack2jack(o, exp, n) 


converts a polynomial expression involving Jack 
polynomials into a linear combination of 
Jack polynomials 


jackidentity 


jackidentity(a, k, m) 


the Jack polynomial evaluated at Xi = 1, i = l..m 


jacobi 


]&coh\{a,K,gi,g 2 , [x]' N') 


the multivariate Jacobi polynomial as a linear 
combination of Jack polynomials 


laguerre 


laguerre(a, n, g, [x],' N') 


the multivariate Laguerre polynomial as a linear 
combination of Jack polynomials 


leg 


leg(K,i,j) 


the leg-length of a partition at the square (i,j) 


lhook 


lhook(a, k) 


the lower hook of a partition 
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Procedure Syntax Description 



m2jack 


m2jack(o!, exp, n) 


converts an expression involving monomials in 
Jack polynomial basis 


m2m 


m2m(exp, n) 


converts an expression involving monomials 
to a linear combination of monomials 


m2p 


m2p(exp) 


converts an expression involving monomials 
to a linear combination of power sum functions 


p2m 


p2m(exp, n) 


converts an expression involving power sum functions 
to a linear combination of monomials 


par 


par(fc) 


produces and lists all partitions of a given integer 


rho 


rho(a, k) 


the p function of a partition 


sfact 


sfact (r, k) 


the shifted factorial (Pochhammer symbol) 


subpar 


subpar (n) 


produces and lists all subpartitions of a given partition 


uhook 


uhook(a, k) 


the upper hook of a partition 



Figure 2: Dependence graph for the procedures of MOPS. 



par arm conjugate sfact rho issubpar m2m 




3.4 Computing Expected Values / Evaluating Integrals 

Let P\(x\, . . . ,x m ) be a symmetric polynomial in m variables, with highest 
order term corresponding to the partition A. To compute the expected value 
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of P\ with respect to one of the distributions dfij Q1UJI . dfj," dj) , or cfyx"^ 
(|17j) . we write 



-PaC 3 -!) • • • j x m) — ^ ] C/j^C^Ii, 



and by applying the linearity of expectation, we obtain 

£[P A Oi,- • • ,x m )) = ^2,c K , a E[C%{x\, . . .,x m )] . 

K 

In the univariate case the Jack polynomials are simply monomials, and 
we have the following (well-known) moments for the Hermite, Laguerre, and 
Jacobi weight functions: 

-±= f x k e- x2 ' 2 dx = (2fc-l)H = (-l)*/ 2 fl- fc (0) , 
F(2 + a + 6) f lV(1 _ l)t(fa = (a + l) t r(q + > + 2) j..» (0) . 



r(o + i)r(6 + 1) y [0) i] v ' r(a + i)r{a + b + k + 2) 

In the above, k > 0. 

A similar triad of formulas is can be established for the multivariate case. 
In the Laguerre and Jacobi cases, the univariate formulas generalize easily: 

f rn — 1 

/ C%{x 1 ,...,x m )dnt(x 1 ,...,x m ) = ( 7 + + l)«C*(/ m ) = L%«(p) , (24) 

</[0,oo) m a 

C%(x 1 ,...,z m )dn%x 1 ,...,x m ) = 7^27 ^ = 4 a ' 9l ' g2 (0)(25) 

[o,i]™ {9i+g2 + i{m - 1) + 2) K 

For a good reference for the first formula, see Forrester and Baker [2] ; the 
second one was obtained by Kadell JH] • 
For the Hermite case, 



/ C2( Xl ,...,x m )diJ% = (-i) k / 2 n2(o), 



(26) 



but to the best of our knowledge, no simpler closed-form formula is known. 
We compute the right hand side as the 0th order coefficient of the polynomial 
H%(x\, . . . , Xm ), using formula ([2*2*]). Note that if k sums to an odd integer, 
the above is trivially 0. 

The procedures expHjacks, expLjacks, and expJjacks compute the 
expected value of an expression (allowing not only for addition, but also for 
multiplication and powers) involving Jack polynomials (all having the same 
parameter a as the distribution). They consist of the two steps described in 
the first paragraph of this section: the first one is reducing the expression 
to a weighted sum of Jack polynomials (if the expression is a weighted sum 
of Jack polynomials, this step is skipped), and the second step is replacing 
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each Jack polynomial with its expected value, using the formulas ()26[) . I|24j). 
and (|25|) . 

The procedures expH, expL, and expj compute the expected value 
of an expression involving monomials (allowing for addition, multiplication, 
and powers), and there are three steps involved: the first is to write the 
expression in monomial basis, the second - to rewrite the result in Jack 
polynomial basis, and the third is to replace the Jack polynomials by their 
expectations, using l(2*H|) . (|2*1|) . and l(23|) . 

Example. Suppose we want to compute the expected value of 

z(a,x 1 ,x 2 ,x 3 ) := J^(x 1 ,X2,x 3 )C^ 1 ^(xi,x 2 ,x 3 ) 

over the 2/a-Hermite distribution. First we have to express z as a linear 
combination of Jack "C" Polynomials. Note that the number of variables, 
as well as a, must be the same in the two terms of z. 
First, we express the two terms in monomial basis: 

Jp A] (x 1 ,x 2 ,x 3 ) = (2 + a) m[ 2 ,i](xi,x 2 ,x 3 ) +6 m[ 1AA ](x 1 ,x 2 ,x 3 ) , 

6a 2 

C\{ AA] (xi,x 2 ,x 3 ) = ( 1 + a ^ 2 + a ) m [i,l,i]( a; i» a; 2» a; 3)- 

Their product thus becomes a linear combination of sums of products 
of two monomials, which are in turn converted in linear combinations of 
monomials. Note that here we use the fact that there are three variables: 

m[2 i i](x 1 ,X2,x 3 ) m[i tljl ] (x 1 ,x 2 ,x 3 ) = m[ 3i2t i](x 1 ,X2,x 3 ) , while 

m[ hlj l](x 1 ,X 2 ,X 3 ) 2 = I7l[ 2i 2,2](si, x 2 ,x 3 ) . 

Putting it all together, in monomial basis, 

. . 6a 2 , 

z{a,x 1 ,x 2 ,x 3 ) = — — rn[ 32 ,i]{xi,X2,x 3 ) + 
1 + a 1 

36a 2 

+ (l + a)(2 + a) m PA2](^i^2,x 3 ). 

All that is left now is to convert from monomial basis back to Jack poly- 
nomial basis. We obtain that 

1 (2 + 3a)(l + 2a) 2 na 

z{a,x 1 ,x 2 ,x 3 ) = — — — r C, 321] {x 1 ,X2,x 3 ) 

lz(J a(l + a) 1 > ' 1 

We are now able to finish the work: 

36(a-l)(a + 3) 



E H [z(a,x 1 ,x 2 ,x 3 )] 



(l + a)(2 + a) 



25 



4 Complexity bounds and running times 



In this section we will analyze the performance of the main algorithms, which 
we divide into four parts: 

1. algorithms that compute polynomials; 

2. algorithms that evaluate integrals; 

3. conversion algorithms; 

4. numerical algorithms. 

Our complexity bounds are upper bounds, but we believe many of them 
to be asymptotically correct. They work well for the numerical evaluation 
of the parameters involved (i.e. a, m, 7, g%, #2); symbolic evaluation of the 
polynomials is considerably slower. We are not aware of the existence of a 
good symbolic performance model for Maple, and hence it would be difficult 
to predict how much slower symbolic evaluation is than numerical evaluation. 
Once parameters are introduced (like m, the number of variables, or a, the 
Jack parameter), the quantities to be computed become rational functions of 
these parameters, of degrees that can go up to the partition size |«|. Storage 
then becomes an issue, hence one would expect that the running times for 
symbolic evaluation would be orders of magnitude slower than for numerical 
evaluation, since the coefficients we deal with must be written and stored 
on "slow" memory (e.g. disk space), and the "transport" time to and from 
"slow" memory greatly increases the overall running time. 

For each algorithm we provide a complexity analysis, and we illustrate the 
performance in practice by providing running times for different tests (both 
numerical and symbolic); then we examine the running times and draw a set 
of conclusions. 

Each time we use N/A for an entry in a running times table, we have done 
so because that particular computation has exhausted the memory available 
to Maple, and hence (regardless of the time it took up to that point) the 
computation was not finished. 

The computer on which we have performed our tests is a Pentium 4 by 
Dell, 1.8 Ghz, 512 MB; the version of Maple used for the tests is Maple 8. 

The last thing worth mentioning is that Maple has an option remember, 
that is it allows for storage and recall of a quantity that was computed 
previously, and that MOPS is taking advantage of that. 

4.1 Algorithms that compute polynomials 

In this category we have the algorithms that evaluate jack, gbinomial, 
hermite, laguerre, and jacobi. We analyze here gbinomial, though it 
is not a polynomial in {x\, . . . ,x m ), because it is the main building block 
for hermite, laguerre, and jacobi, and its complexity determines their 
computational complexity. 

Throughout this section, we will follow the notations given in Table [7| 
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number of partitions of k dominated, by k 




number of partitions of the number k (each partition of 


k is dominated by [k]) 




number of subpartitions of k 




number of subpartitions of k which are superpartitions 




for a (this implies a is a subpartition of k) 




number of subpartitions of k which sum to a number with 




the same parity with k 



Table 7: Notations to be used throughout Section 



To make estimates, we have used Ramanujan's formula: 

1 

4ky/S 



P [k] ~ -L= , (27) 



and the inequalities 

A K < U K < P K < P [k] and U Kj(T <P [k] (28) 

for asymptotical estimates. 

1. jack. The algorithm uses recurrence (jSJ), together with the 'boundary 
conditions' c K) \ = if k )£_ A in dominance ordering, and c K)K = c ° K fc ^ ■ 

The length of the recurrence is at most 0(/ci( fc ^' 1 )), with k\ being the 
first entry in the partition, and the algorithm will check each of the 
possible partitions \i (at most k\ i^^) ) to see if they are dominated by 
k and dominating A (this involves I additions and I comparisons). The 
rest of the computation has complexity 0(k). 
Thus the complexity of the algorithm is 0{k\k ?J P K ). 
Using the inequalities (|28|) . the best asymptotical upper bound we 
can get is for the complexity of computing a Jack polynomial is thus 
0(k 3 e 7T ^ 2k / 3 ), which is super-polynomial. 

Below we illustrate the running times for both numerical and symbolic 
computations. For numerical computations, we have chosen to make 
a = 1, so that the Jack polynomials are the Schur functions. Note that 
we do not test the partition [k] ; for that particular partition we have a 
closed- form formula for the Jack polynomial, due to Stanley [HJ , which 
has complexity 0(kP k ) ~ O^v 7 ^). 
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T^mrninp" tiinp rv = 1 

A ^ l L JL_L J. J. 3. J. J. £^ U XX±±Vj ^ 


Rnnninp^ timp rv svmholir 


a t i o 


1 5 


K = 




2.48 


4 54 


1 83 




K — — 


8 71 
L°> 'J 


1.79 


3.17 


1.77 




K" 

rv — 


ro o o o ol 
[O, O, O, O, OJ 


0.39 


0.50 


1.28 


2D 


K = 


[19,1] 


1 6 97 


30 45 


1 79 




K = 


[10, 10] 


11.53 


20.32 


1.76 




K = 


[4,4,4,4,4] 


2.91 


4.02 


1.38 


25 


K = 


[24,1] 


93.42 


189.66 


2.03 




K = 


[9,8,8] 


46.85 


79.85 


1.70 




K = 


[5,5,5,5,5] 


16.08 


24.18 


1.50 


30 


K = 


[29,1] 


634.32 


1819.65 


2.86 




K = 


[10, 10, 10] 


214.10 


418.19 


1.95 




K = 


[6,6,6,6,6] 


73.54 


113.55 


1.54 



Table 8: Running times (in seconds) for the Jack polynomial computation. 

Remark 4.1. Note that the ratio of the running times increases when 
the partition size increases. At k = 30, the number of partitions is 
5604, and each of the monomial coefficients is a rational function of a. 
Issues like storage and memory access become important, and influence 
negatively the running times. Another important factor is that in order 
to make things easier to store and access, not to mention easier to read 
and interpret, we use the procedures "simplify" and "factor", which are 
relatively costly. 

Extrapolation. Since the speed/memory of a top-of-the-line computer 
seems to go up by a factor of 10 3 every 10 years, one can predict that 
within a decade, using MOPS, computing JZ g ^ will take about 30 
minutes. 

2. gbinomial. We use together with the boundary conditions listed 
in Section \2. 51 and with the contiguous binomial coefficient formula Q. 
From (JHJ), it follows that computing a single contiguous binomial coef- 
ficient has complexity 0(k), and one needs to compute no more than I 
such coefficients per subpartition a of k which is a superpartition of a. 
Thus one immediately obtains the bound 0{klU Ka ) for the complexity 
of computing This is smaller than 0(k 2 U K ^ rpi). 

Note that by computing one also obtains for each a C fj, c k. 
So we have chosen for our tests to compute (rjjn) f° r different k, as this 
yields all the binomial coefficients having k as top partition (except (Z) , 
but that requires only an additional complexity O(kl)). 
By using the inequalities (|28[). we obtain an asymptotical upper bound 
of 0(ke n ^ 2k ^) for computing all the generalized binomial coefficients 
corresponding to partitions of k. 
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k 


K 


Running time, 


Running time, 








a = 1 






1 5 


[6 A 9 9 11 

[U, 1, z, Z, 1J 


99 


1 19 


1 

IOC/ 




3 3 3 3 31 


0.05 


0.18 


56 




10 51 


0.03 


0.15 


51 




[64 3 9 9 1 1 11 

[U, 4:, O, Z, Z, 1, 1, Ij 


1 01 


6 68 


41 8 

IIO 




4, 4, 4, 4, 41 


0.17 


0.6 


126 




[12,8] 


0.07 


0.28 


81 


25 


[7,5,4,3,2,2,1,1] 


3.41 


23.37 


1077 




[5,5,5,5,5] 


0.41 


1.67 


252 




[16,9] 


0.15 


0.62 


125 


30 


[8,6,4,3,2,2,1,1,1,1,1] 


11.87 


89.61 


2619 




[6,6,6,6,6] 


0.91 


3.95 


462 




[20, 10] 


0.24 


1.20 


176 



Table 9: Running times (in seconds) for the generalized binomial coefficient computa- 
tion. 

Remark 4.2. Once again, size and length of the partition increase the 
symbolic running times; however, note that the running times are rela- 
tively small, even for partitions of 30. We believe that the generalized 
binomial coefficients are rational functions of a which can always be fac- 
tored in small-degree factors, so that they are easy to store and operate 
with. 

3. jacobi. 

To compute the Jacobi polynomials, we use the format of Section f2. 7. 21 
and recurrence (|14jl. One can easily see that at each step, one needs to 
compute at most I contiguous binomial coefficients, each of which has 
complexity 0(k); in addition, one needs to compute another at most I 
binomial coefficients; each of these takes only 0(1), as the contiguous 
coefficients needed have already been computed at the previous step. 
Thus the total complexity is 0{kl) (since I < k) at each step, for a total 
ofO(M£7-« >[ia] ). 

Hence computing numerically the Jacobi polynomials is comparable 
to computing the generalized binomial coefficients (m^i) j however, the 
constant for the Jacobi polynomial complexity is considerably larger 
(our best guess sets it around 8). 

The best asymptotical upper bound we can obtain using the inequalities 
(t28l) is thus once again 0(ke 1 "^^). 

The Jacobi parameters we chose for each of the computations below are 
and 1. 

Remark 4.3. While the running times for numerical evaluation are 
reasonable, they explode when a symbolic parameter is introduced. The 
coefficients of the polynomial are rational functions of that parameter 
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k 




Running time, 


Running time, 


Running time, 








£J = 1, TYl = I 


777, symbolic 


at, m symbolic 




10 


[4 2 2 1 11 


0.27 


0.74 


22.12 


42 




4 3 31 


0.11 


0.35 


1.88 


30 




F7 31 


0.10 


0.30 


1.57 


26 


15 


[fi 4 2 2 11 


1.05 


11.08 


N/A 


139 




3,3,3,3,3 


0.39 


0.87 


63.07 


56 




[10,5] 


0.19 


1.01 


27.98 


51 


20 


[6,4,3,2,2,1,1,1] 


5.94 


N/A 


N/A 


418 




[4,4,4,4,4] 


0.63 


8.24 


N/A 


126 




[12,8] 


0.26 


3.51 


N/A 


81 


25 


[7,5,4,3,2,2,1,1] 


18.61 


N/A 


N/A 


1077 




[5,5,5,5,5] 


1.23 


N/A 


N/A 


252 




[16,9] 


0.45 


N/A 


N/A 


125 



Table 10: Running times (in seconds) for the Jacobi polynomial computation. 



or combination of parameters, of order up to k{k — l)/2. We recall that 
there are U Kj rpi of them, a potentially superpolynomial number, which 
explains the tremendous increase in the running time. 

4. laguerre. 

We use the format given in Section 12.7.41 it is easily established that 
the complexity of computing the Laguerre polynomial is dominated by 
the cost of computing the binomial coefficients, that is 0(klU K ^ h?]), 
and once again the best asymptotical upper bound we can obtain using 
the inequalities (|28|) is thus once again 

The Laguerre parameter we chose for each of the computations below 
is 1. 

Remark 4.4. For the Laguerre polynomials, even in the all-symbolic 
case, the computation is very easy, and the storage required is relatively 
small. This explains why it is possible to obtain them without much 
effort, in any one of the cases. 

5. hermite. 

We use the format given in Section 12.7.61 and recurrence (|2*T|) . We 
only do work for those coefficients that correspond to subpartitions a 
of k such that |er| = k (mod2). There are A K of them. For each, we 
compute at most (2) contiguous coefficients, each computed with 0(k) 
complexity. The complexity of the rest of the computation is 0{k). 
Hence the total complexity is 0(kl 2 A K ). 

Remark 4.5. A K = 0(U K ); A K ~ U K /2. 

Hence one asymptotical upper bound that we can obtain for the com- 
plexity of computing a Hermite polynomial is 0(k 2 e 7T ^ 2k ^ 3 ). 
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k 


K 


Running time, 


Running time, 


Running time, 


^ K 






r\i — 1 fji — / 

(_X X , lib li 


lit o v iiiuuiHj 


r\ fn t;VTTll~inllP 




i n 


\A 9 9 1 11 

[% Z, Z, 1, 1J 


1 9 


9^ 


n 54 


L ±Ai 




4 3 31 


0.07 


0.14 


0.31 


30 




7 31 


0.07 


0.10 


0.28 


26 


i 5 


\(\ 4 9 9 11 


49 


89 


9 95 


1 S9 




3 3 3 3 3 


0.18 


0.27 


0.84 


56 




[10,5] 


0.11 


0.22 


0.81 


51 


20 


[6,4,3,2,2,1,1,1] 


2.26 


3.37 


16.08 


418 




[4,4,4,4,4] 


0.44 


0.69 


2.74 


126 




[12,8] 


0.20 


0.37 


1.79 


81 


25 


[7,5,4,3,2,2,1,1] 


7.23 


11.06 


67.92 


1077 




[5,5,5,5,5] 


0.96 


1.53 


8.06 


252 




[16,9] 


0.32 


0.69 


4.21 


125 


Table 11: Running times (in seconds) for the Laguerre polynomial computation. 


k 


K 


Running time, 


Running time, 


Running time, 


A K 






a = 1, m = I 


m symbolic 


a, m symbolic 




10 


[4,2,2,1,1] 


0.21 


0.24 


0.75 


22 




[4,3,3] 


0.09 


0.11 


0.33 


16 




[7,3] 


0.05 


0.06 


0.24 


14 


15 


[6,4,2,2,1] 


0.41 


2.83 


42.92 


88 




[3,3,3,3,3] 


0.13 


0.17 


1.83 


38 




[10,5] 


0.10 


0.12 


1.10 


30 


20 


[6,4,3,2,2,1,1,1] 


1.93 


2.39 


N/A 


211 




[4,4,4,4,4] 


0.35 


0.51 


N/A 


66 




[12,8] 


0.18 


0.25 


13.49 


43 


25 


[7,5,4,3,2,2,1,1] 


6.23 


7.53 


N/A 


1077 




[5,5,5,5,5] 


0.90 


1.20 


N/A 


252 




[16,9] 


0.29 


0.50 


106.56 


125 



Table 12: Running times (in seconds) for the Hermite polynomial computation. 

Remark 4.6. Note that when m is parametrized, but a = 1, the 
computation is almost as fast as in the all-numerical case. That hap- 
pens because the dependence on m is very simple, and it only involves 
Pochhammer symbols, which do not get expanded (so that the storage 
required is minimal). However, the dependence on a is more compli- 
cated, and the rational functions obtained as coefficients are complex 
and hard to store. Hence the running time for the all-symbolic compu- 
tation increases dramatically. 
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4.2 Conversion algorithms 

There are five conversion algorithms, m2jack, jack2jack, m2m, p2m, and 
m2p. 

1. m2jack. This algorithm computes and then inverts the change of ba- 
sis matrix from monomials to Jack polynomials, taking advantage of 
the fact that the matrix is upper triangular. At each turn, the algo- 
rithm extracts the highest order monomial remaining, computes the 
coefficient of the corresponding Jack polynomial, and then extracts the 
monomial expansion of the Jack polynomial from the current monomial 
expression. 

Let k be the highest-order monomial present in the initial expression, 

and let us assume that the expression is homogeneous. 

Then the complexity of the computation is dominated by the complexity 

of computing the Jack polynomial expansion in terms of monomials for 

all partitions of k smaller in lexicographical ordering than k. 

It follows that an upper bound on the complexity is given by 0(D K k A D K ) = 

0{k 2 e 2 *^^). 



The performance in practice is exemplified below. 



Partition sum 


Partition 


Runtime 


Runtime 


Ratio 






(a = 2) 


symbolic a 


of the two 


fc = 6 


«=[6] 


0.14 


0.45 


0.31 


k = 7 


K=[7] 


0.29 


1.04 


0.27 


k = 8 


K=[8] 


0.63 


2.91 


0.21 


k = 9 


K=[9] 


1.21 


7.49 


0.16 


k = 10 


K = [10] 


2.62 


20.25 


0.12 


k = 11 


« = [11] 


4.77 


54.75 


0.08 


k = 12 


K = [12] 


8.82 


186.09 


0.04 


fc = 15 


K = [15] 


52.65 


7177.02 


< 0.01 



2. m2m. The algorithm takes an expression involving products of mono- 
mial functions and writes it in monomial basis by deciding which par- 
titions appear in the expansion and by counting the number of times 
they appear. Hence this algorithm is an alternative to adding the m 
basis as the dual of h in the SF package, and using the torn procedure 
afterward (though the torn procedure is more general than this). 
We have tested m2m against torn, and we have found that on partitions 
where the ratio sum-to-length is high, m2m performs much better, 
while on partitions with the ration sum-to-length is small, the tables 
are reversed. Hence we recommend to the user who wants to use our 
library, but might be working with partitions of the latter case, to also 
obtain and install SF and use it for computations. 

Below is a performance comparison. The number of variables n used 
in m2m, each time, was the sum of the partition lengths (which is the 
smallest number of variables that requires obtaining all terms). 
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Input 


n 


Runtime m2m 


Runtime torn 


Ratio 


m[5,2,l,l] -m[4] 


5 


0.07 


0.28 


0.25 


m[3,2,l] -m[5] 


4 


0.03 


0.10 


0.30 


m[5,3,2] -m[4,3] ■ m[2] 


6 


3.52 


35.27 


0.10 


m[7,3, 1] -m[4,2] 


5 


0.30 


9.89 


0.03 


m[4,3,2] -m[6,4] 


5 


0.29 


35.72 


< 0.01 


m[2,2,l] -m[l,l] 


5 


0.05 


0.03 


1.66 


m[3,l] -m[2,2,l] 


5 


0.12 


0.05 


2.40 


m[2, 1,1,1] -m[2,l] 


6 


0.22 


0.04 


5.50 


m[3,2, 1] -m[2, 1,1,1] 


7 


2.95 


0.10 


29.5 


m[3, 1,1] -m[2,l] -m[l,l,l] 


8 


13.55 


0.13 


104.23 



3. p2m. The algorithm expands a product of simple power sum functions 
into monomial basis. This is an alternative to adding the m basis as the 
dual of h in the SF package, and calling the torn procedure with power 
sum functions as inputs. As was the case with m2m, our algorithm 
performs much better on partitions with high sum-to-length ratio, and 
torn performs better on partitions with low sum-to-length ratio, as can 
be clearly seen from the performance comparison below. 



Input 


Runtime p2m 


Runtime torn 


Ratio 


PA ■ P>\ ■ P2 


0.03 


0.14 


0.21 


P8-P5- P{ 


0.20 


5.46 


0.04 


Pi -Pa- Ps 


0.01 


0.46 


0.02 


P\-P\- P2 


0.04 


5.40 


< 0.01 


PA-P2- Pi 


0.12 


0.10 


1.20 


P5-P2- Pi 


0.96 


0.17 


5.64 


P3-P2- Pi 


1.97 


0.04 


49.25 


A A 


16.27 


0.15 


108.46 



4. m2p. The algorithm converts an expression of monomials into power 
sum functions; it is an alternative to the top option in the SF package. 
As before, for high sum-to-length ratio, our algorithm performs better, 
whereas the reverse is true for low sum-to-length ratio. It is perhaps 
worth noting that for this case, the ratio sum-to-length has to be be- 
tween 1 and 2 for a significant outperformance of our m2p by top to 
occur. This can be seen in the performance examples below. 



Input 


Runtime m2p 


Runtime top 


Ratio 


m[l, 1,1, 1,1, 1,1,1] 


3.61 


0.04 


90.25 


m[3,2,2,2,l,l] -m[2, 1,1,1] 


2.19 


0.11 


19.91 


m[2,2, 1,1,1] -m[2,l] 


0.120 


0.03 


4.00 


m[l, 1,1,1] ■ m[l, 1, 1] 


0.03 


0.02 


1.50 


m[4,3,2,2] 


0.06 


0.18 


0.33 


m[10, 1] 


0.01 


0.10 


0.10 


m[5,4,3] -m[3] 2 


0.02 


0.23 


0.08 


m[5,4,3,3] 


0.08 


3.21 


0.02 


m[3,3,2, 1] -m[5,2] 


0.07 


5.57 


0.01 
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5. jack2jack. This algorithm takes an expression in Jack polynomials, 
and turns it into a linear combination of Jack polynomials, by tak- 
ing each multiplicative term, expanding it in monomial basis, then us- 
ing m2m to get rid of the resulting multiplicative factors, and finally, 
m2jack to convert the linear combination of monomials back in Jack 
polynomial basis. 



Input 


n 


Runtime (symbolic) 


Runtime (a = 1) 


J[2,l] -C[l,l,l] 


3 


0.04 


0.03 


J[3,2,l]-C[3] 


5 


3.53 


0.45 


C[2f • J[4,2] 


5 


11.44 


1.17 


C[2Y- J[3,l] 


8 


29.09 


3.35 


P[3,2] • J[2,l,l] 


7 


15.00 


2.10 


P[2,1].J[4,2]-C[2] 


5 


28.95 


2.07 



4.3 Algorithms that evaluate integrals 

Here we have expHjacks, expLjacks, expJjacks, expH, expL, and expJ. 

These algorithms depend on the length and complexity of the input. Let 
P be the polynomial one wishes to analyze; one must first convert P to a lin- 
ear combination of Jack polynomials, and then replace each Jack polynomial 
with its expectation. 

Case 1. Suppose P is in monomial format, as an expression which involves 
sums and products of monomials. First we convert P to a linear combination 
of monomials using m2m, and then we convert that linear combination of 
monomials to a linear combination of Jack polynomials using m2jack. 

For any term of the form uijimj! . . . m\ P , with A , A 2 , . . . , X p not neces- 
sarily distinct partitions, when we expand it in monomial basis, the largest 
possible number of terms is D^, where fi is the partition which results 
from the superposition of A 1 , A 2 , . . . , A p , i.e. fi\ = \\ + A 2 + . . . + A p , 
f-2 = A2 + A 2 . + . . . + Ag, etc.. Let u = 

After the expansion in monomial basis, applying m2jack on the resulting 
expression has complexity 0(u 4 D^) = 0(n 2 e 37r v /2 /3v / "), 

Remark 4.7. This however is a very relaxed upper bound, and if we start off 
with P being a sum of a few (n) monomials, the call to m2m is not executed, 
and the complexity of the call to m2jack is 0(nn 4 L> 2 ) = 0(nu 2 e 37r v /2 7V"). 

As explained in Section 13.41 the first step is common to expH, expL, 
and expJ. The second step is different and its complexity is much higher 
for expH than for expL or expJ. However, as we can see from the running 
times in the table below, the calls to m2m and m2jack (made in the first 
step) are much more expensive than the substitutions, and so the overall 
running times are comparable. 

In these examples, we consider a symbolic parameter a, a symbolic num- 
ber of variables n, 7 = 1, and g\ = gi = 1. 



34 



Input 


Runtime expH 


Runtime expL 


Runtime expJ 


m[6] 


0.94 


0.70 


0.80 


m[3,3,2] 


1.98 


0.85 


0.96 


m[5,2,l] 


5.69 


3.20 


4.20 


m[3, 1,1,1] -m[2] 


4.23 


1.84 


2.59 


m[4, 1] -m[l,l,l] 


3.94 


2.18 


3.58 


m[5, 1] • m[2] 


8.86 


6.04 


9.82 


m[3] 2 • m[2] 


8.80 


7.00 


13.04 


m[4,2] -m[3, 1] 


39.85 


35.71 


68.82 



Case 2. Suppose P is in Jack polynomial format; then we use jack2jack 
to write the it as a linear combination of Jack polynomials, and finally we 
replace each Jack term by its expected value. The first step, as before, is 
common to all three procedures (expHjacks, expLjacks, expJjacks). 

While in the case of expHjacks the complexity of computing the ex- 
pectation is 0(-u 4 e 27r \^7V"), in the cases of expLjacks and expJjacks 
the same complexity is only 0(u). This explains the significant differences 
recorded in the first three rows of the table. It is also worth noting that in 
the case of an odd u, the time it takes to compute the expected value of 
a Jack polynomial with Hermite weight is 0, as the value of the output is 
known in advance to be 0. 

The complexity of expressing a product of Jack polynomials in Jack poly- 
nomial basis is much higher than the computation of a single Jack polyno- 
mial expected value. This explains why, in the last few rows of the table, 
the entries are no longer so different in magnitude. 



In the examples below, we considered a symbolic parameter a, a symbolic 
number of variables n, 7 = 1, and g\ = gi = 1. 



Input 


Runtime expHjacks 


Runtime expLjacks 


Runtime expJjacks 


C[4,3,2,l] 


0.30 


0.03 


0.03 


C[6,6,2] 


1.06 


0.04 


0.04 


C[7,5,3,l] 


4.70 


0.04 


0.05 


C[10,3,2, 1] 


4.47 


0.05 


0.05 


C[3, 1]-P[2,2] 


14.75 


12.75 


12.93 


C[4,2].J[1,1] 


31.86 


29.05 


30.11 


J[2, 1, 1] • J[4] 


76.62 


81.93 


80.14 


C[2,2].J[4] 


53.79 


54.30 


55.07 



4.4 Numerical algorithms 

Some of the symbolic/numerical evaluation routines analyzed in the pre- 
vious sections include options for polynomial evaluation on numerical val- 
ues of the x variables. The routines that compute the polynomials Jack, 
Hermite, Laguerre, and Jacobi have options that allow for numerical val- 
ues of the x variables. This makes it possible to compute quantities like 
Cjg 2 j (2.53, —1.09, 7.33); this feature can be used for graphics (when one needs 
to plot some statistic of a random matrix, as we demonstrate in the next 
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section). 

The algorithms we have used to implement these options have been de- 
veloped and analyzed by Koev and Demmel [0] for the Jack polynomials; to 
evaluate the other polynomials, we use the regular expansion in terms of Jack 
polynomials, then substitute the numerical values for each Jack polynomial. 

5 Applications 

We have written this library for the user who would like to do statistical 
computations, form or test conjectures, and explore identities. The great 
benefit is that all computations can be done symbolically, keeping a as a 
parameter; the downside of symbolic computations, as we have mentioned 
before, is that the storage space required is very large, and computations are 
consequently slowed down. Our experience, however, was that on a good, 
but not top-of-the-line machine (see specifications in Section 0J, we have 
been able to increase the size of the partition enough in order to make and 
then satisfactorily test conjectures. 

Below are some examples of computations that we imagine are of the type 
a researcher might want to use in forming conjectures, or of the type that 
might be useful in practice. 

Some of the applications, like the computation of the moments of the 
trace, can be done with symbolic a and n (number of variables); others, 
like the computation of the moments of the determinant, need an actual 
value for n, but allow for symbolic a computations; yet others, like the level 
density computation, need all numerical parameters. For each computation, 
we have tried to indicate upper bounds for the size of the necessary numerical 
parameters. 

1. Moments of the determinant. One of the many interesting problems in 
random matrix theory is computing the moments of the determinant of 
a square random matrix. If the eigenvalues are chosen to have the 2/a- 
Hermite distribution (given by the weight function the problem 
of computing the determinant is non-trivial. Closed form answers are 
known for the cases a = 1/2,1, and 2 (see [Q, 0, |2Zj); however, 
the general a case does not have an explicit answer (except for some 
particular situations like in chapter 8]. 

Since the kth moment of the determinant's distribution is given as the 
integral of m^m] (xi, . . . , x m ) = C^ m] (xi, . . . , a? m )/C[ fc m] (I m ) over the 
corresponding 2/a-Hermite distribution, MOPS can be used in evalu- 
ating it for specific values of k and m. 

For example, for k = 2 and m = 5, the answer can be obtained by 
typing in 

> factor (expHjacks (a, C[2,2,2,2,2], 5)/jackidentity(a, [2,2,2,2,2], 5)); 
and the output is 

a 4 + 10 a 3 + 45 a 2 + 80 a + 89 

> a* 
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The duality principle between a and 1/a proved in 7, Section 8.5.2] 
linking the expected value of the kth. power of the determinant of a 
n x n matrix to the expected value of the nth power of a k x k matrix 
is illustrated below: 

> factor (expHjacks(l /a, C[5,5], 2)/jackidentity(l/a, [5,5], 2)); 

with output 

> - a (a 4 + 10 a 3 + 45 a 2 + 80 a + 89) 

Remark 5.1. In practice, we have observed that computations with 
a symbolic and k ■ m < 22 can be performed relatively fast (under 2 
minutes on the computer with specs given in the beginning of Section 
0); for k ■ m > 22 and a symbolic, the amount of memory available 
begins to play an important role. For actual values of a (for example, 
a = 1 ), the computation for k = 10 and m = 5 took under 40 seconds. 

2. Expectations of powers of the trace. Consider the problem of computing 
the expected value of the 6th power the trace of a Hermite (Gaussian) 
ensemble (here n is an arbitrary integer). This amounts to making a 
call to expH, simplifying, and expanding the answer in Taylor series 
for a clear format. In short, a one-line command: 

> taylor(simplify(expH(a, m [6], n)), n); 
with answer 

15 a 3 - 32 a 2 + 32 a - 15 -54 a + 32 a 2 + 32 2 22 a - 22 3 5 4 

a 3 a 3 a 3 a 3 



Remark 5.2. This computation emphasizes best the power of MOPs. 
It is very quick (took 0.8 seconds on the test machine (see specifications 
in Section^ and it allows for both a and n symbolic. The same com- 
putation for the 12th power of the trace with a and n symbolic took less 
than 8 minutes. 

Integrals of powers of the trace are related to Catalan numbers and 
maps on surfaces of various genuses, and are of interest to (algebraic) 
combinatorialists ( ) . 

3. Smallest eigenvalue distributions. One of the quantities of interest in the 
study of Wishart matrices 6 is the distribution of the smallest eigenvalue. 
There is an extensive literature on the subject, starting with the work of 
James |16j and Constantine [3]. More recent references are Silverstein 
30 and Edelman pj. In |7j, we find a closed- form answer for general a 
and integer values of 7, in terms of a hypergeometric 2-^0 function (see 
also (HU). 

We wrote a small script (presented below) implementing the formula, 
and used it to compute the exact distribution of the smallest eigenvalue 
of a Wishart matrix for a = 1 (the complex case) for n = 3, m = 6, and 
n = 2, m = 10, which we plotted in MATLAB. We have also used a 

6 The joint eigenvalue distribution of Wishart matrices is given by the Laguerre weight fi^' 1 with 
a = 1 (complex case) or a = 2 (real case). 
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Monte Carlo simulation to plot in MATLAB histograms of the smallest 
eigenvalue of matrices from the corresponding Wishart ensemble, for 
comparison (see Figure 0J- For the histograms, we have chosen in each 
case 30, 000 samples from the corresponding Wishart ensemble. 

smalleig:=proc(n,k,x) local r,t,i, y,inte; 
if (n>l) then r: = [-2/x]; 
end if; 

for i from 2 to (n-1) do 

r: = [op(r),-2/x]; 
end do; 

t:=x~ ((k-n)*n) * exp(-x*n/2) * ghypergeom(l, [n-k, n+l],[],r,'m'); 
return simplify (t); 
end proc; 

scaledsmalleig:=proc(n,k,x) local inte, yy, z; 

yy :=z->smalleig(n,k,z); 

inte := integrate(yy(z), z=0.. infinity); 

return(smalleig(n,k,x) /inte); 
end proc; 



zz:=scaledsmalleig(3,6, x); 
plot(zz, x=0..10); 



n=3, m=6, P = 2 n=2, m=1 0, p = 2 




Figure 3: Histograms of the smallest eigenvalue distribution for the complex Wishart 
ensembles of size (3,6) and (2,10) (a = 1), together with the exact distributions as 
given by p3|) . 

4. Level densities. Level density formulas are well-known in terms of or- 
thogonal polynomials for a = 1/2,1,2. Forrester and Baker [3] have 
computed these densities in terms of a multivariate Hermite polynomial 
for P = 2/ a an even integer (i.e. a is the inverse of an integer). We 
have found an equivalent formulation for the level density of the n x n 
Hermite ensemble for which a is the inverse of an integer (equivalently, 
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j3 = 2/ a is an even integer). This formula is presented below: 

= ^(- 1 ) n/a ^f e - x2/2 ^/«)"-](^)' 

where the partition [(2/a) n_1 ] is the partition that consists of 2/a re- 
peated n — 1 times. 

To compute the Hermite polynomial, we used the formula (|22[) . and 
in order to get all the eigenvalues roughly in [—1, 1] we scale both the 
variable and the density function by \/2n(3 = \J Anja (see the script). 
We have used the script below to produce Figure 01 which is an exact 
plot of the level densities for n = 4, and (5 = 2, 4, 6, 8, 10 (equivalently, 
a = 1,1/2,1/3,1/4,1/5). 

leveldens:=proc(a,k::list, n, x) option remember; 
local s,u,ut,ul,ks,ss,j,i,sp, result, t,tl,r,jp,ull,c,bbb; 
if(not('MOPS/parvalid'(k))) then return; 
end if; 

result:=0; ks:=sum(k[i] ,i=l..nops(k)); sp: = 'MOPS/subPar'(k); 

we compute the Hermite polynomial evaluated at xl n , using formula 

for s in sp do 
ss:=0; c:=0; ss:=sum(s[i],i=l..nops(s)); 
if not((ss mod 2) = (ks mod 2)) then next; 
end if; 

for j from ss to (ks+ss)/2 do 
jp:='MOPS/Par'(j); 

ull: = (convert(jp,set) intersect convert(sp,set)); 
ul:=Q; 

for ut in ull do 

if 'MOPS/subPar?'(s,ut) then ul: = [op(ul),ut]; 

end if; 
end do; 
t:=0; 

for u in ul do 

tl: = 'MOPS/GSFact'(a,r+(n+a-l)/a,k)/'MOPS/GSFact'(a,r+(n+a- 
l)/a,u); 

t:=t+'MOPS/GBC'(a,k,u)*'MOPS/GBC'(a,u,s)*coeff(tl,r,(ks+ss)/2- 

j); 

end do; 
c:=c+t*(-l)-j; 
end do; 

bbb:=factor(c*(-l) " (ss/2)*x" (ss)); 
result —result +bbb; 
end do; 

result— result*(-l)~(ks)*(-l)~(ks/2) * exp(-x"2/2) * l/sqrt(2*Pi); 
result— result * factor (GAMMA(l + l/a) /GAMMA(l+m/a)); 
end proc; 

9^9^ we scale both the variable and the density function by 
z— (x,b)->sqrt(2*4*b)*leveldens(2/b, [b,b,b], 4, x*sqrt(2*4*b)); 

plot(z(x,2), z(x,4), z(x,6), z(x,8), z(x,10), x=-1.2..1.2, y=-.1..1.4); 
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Figure 4: Level densities for n = 4, a = 
decreases. 



1,1/2,1/3,1/4,1/5; "bumps" increase as a 



1 .4 




-1 



-0.5 



0.5 



x 



For illustration purposes, here is the exact (scaled) density for a = 1/4 
and n = 5, plotted above: 



50685458503680000^77 



(2814749767106560000000000000000 x 32 - 2814749767106560000000000000000 x 30 + 
1720515795143884800000000000000 x 2S - 696386684568207360000000000000 x 2e + 
194340604354756608000000000000 x 2A - 36625240845346406400000000000 x 22 + 
4740055701777285120000000000 x 20 - 658121972672102400000000000 x 18 + 
162266873453346816000000000 x 16 - 31084533121233715200000000 x 14 + 
2673909486122434560000000 x 12 - 136819200341311488000000 x 10 + 
29341248756019200000000 a; 8 - 1130060455927603200000 x G + 
67489799891754240000 x 4 - 2060099901411552000 x 2 + 32632929952848225). 

5. Conjectures. We present here a conjecture that we formulated with the 
help of MOPs. This conjecture was proved later by Richard Stanley. 

Conjecture 5.3. Let k be an integer, a a positive real, and consider 
the representation of the monomial function 



m 



[k] - •/*>" C A • 
X\-k 



Then for all A 




where n(A) is an integer which does not depend on a. 
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6 Copyleft 



Copyleft 2004 Ioana Dumitriu, Alan Edelman, and Gene Shuman. 

Permission is granted to anyone to use, modify, and redistribute MOPs 
freely, subject to the following: 

• We make no guarantees that the software is free of defects. 

• We accept no responsibilities for the consequences of using this software. 

• All explicit use of this library must be explicitly represented. 

• No form of this software may be included or redistributed in a library 
to be sold for profit without our consent. 
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